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By employing Hopf's functional method, we find the exact characteristic functional for a simple 
nonlinear dynamical system introduced by Orszag. Steady-state equal-time statistics thus obtained 
are compared to direct numerical simulation. The solution is both non-trivial and strongly non- 
Gaussian. 



in 
o 
o 

(N 
C 

Hi 



(N 
> 

O 

\o 
o 
in 
o 



> 

x 



PACS numbers: 05.45.-a, 47.52. +j, 05.45.Ac, 05.45.Pq 



I. INTRODUCTION 

A nonlinear dynamical system can be completely deter- 
ministic and its solution unique, yet two trajectories that 
begin close to one another may diverge significantly in a 
finite time. Such sensitive dependence on initial condi- 
tions sets a fundamental limit on predictive accuracy, as 
these systems forget their initial conditions after a short 
time. Probabilistic descriptions, by contrast, avoid the 
details of time-evolution, and instead answer meaningful 
statistical questions. As a canonical example, experimen- 
tal measurements of turbulent velocity fluctuations show 
disordered and unpredictable behavior yet reproducible 
statistical properties. The complexity of the flow con- 
trasts with the smoothness of averaged quantities like 
the energy spectrum, demonstrating that statistical de- 
scriptions may be both economical and insightful. 

To accumulate statistics, numerical simulations may 
implement ensemble averages over initial condition 
and/or long-time integration, yet the computational 
work required to carry out such calculations can be pro- 
hibitive in the case of high-dimensional systems. Purely 
statistical approaches can be quite helpful for such prob- 
lems. Furthermore, sub-grid physics is often modeled 
statistically; fewer inconsistencies will arise if all of the 
dynamical variables arc treated statistically from the out- 
set. Finally, from a theoretical perspective it can be more 
enlightening to work directly with direct statistical ap- 
proaches [1]. 

In this paper we apply a method pioneered by Hopf [2] 
to determine the exact characteristic functional of a sim- 
ple dynamical system introduced by Orszag [3]. Knowl- 
edge of the characteristic functional yields all of the 
equal-time moments and offers some advantage over the 
more common statistical description in terms of the prob- 
ability density function (PDF). In particular, moment 
and cumulant expansions arise naturally in the Hopf for- 
malism [4]. Such expansions, however, suffer from the 
usual closure problem and break down when the statis- 
tics are strongly non-Gaussian. The exact solution to 



Hopf's equation is both non-trivial and, for a finite num- 
ber of dimensions, non-Gaussian. Exact solutions are 
invaluable as they both yield physical insight and pro- 
vide important benchmarks against which approximate 
methods can be tested. 

The outline of the paper is as follows. In Section II we 
introduce the finite-dimensional Orszag-McLaughlin dy- 
namical system. We briefly review the Hopf functional 
approach in Section III. The exact characteristic func- 
tional of the Orszag-McLaughlin system is presented in 
Section IV. In Section V we consider the limit of large 
dimension and show that the statistics become Gaussian 
in this limit. Comparison of the solution against direct 
numerical simulation is made in Section VI. Finally in 
Section VII we highlight a few points. 



II. ORSZAG-MCLAUGHLIN DYNAMICAL 
SYSTEM 

In Orszag's lecture notes [3] and in a subsequent pa- 
per with McLaughlin [5] , a finite- mode system was intro- 
duced as a toy model of inviscid flow: 



dxi 
~dt 



Xi+lXi+2 + Xi-lXi-2 - 2Xi +1 Xi-!, 



(1) 



with 



= 1,.,.,N and periodic boundary condition 
%i+N = The dimension N can be either even or 
odd; we focus on odd N here. Like Euler's equation, 
the equations of motion (EOM) contain only quadratic 
terms. Though the EOM are simple, they generate com- 
plicated trajectories as shown in Fig. 1. 

Despite the fact that the Orszag-McLaughlin dynam- 
ical system is non-Hamiltonian (trivially so in an odd 
number of dimensions), it conserves an energy- like quan- 
tity in the sense that 



E = 



1 N 



(2) 



is a constant of the motion. Trajectories lie on a hyper- 
sphere of radius: 
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FIG. 1: (a) A typical trajectory of the Orszag-McLaughlin 
dynamical system, for the case N = 5, projected onto the x\- 
X2 plane. Trajectories appear random and lie on the surface 
of a 5-dimensional hypersphere. (b) The probability density 
function projected onto the Xi-x<2 plane, as determined by 
long-time integration. 



^ -g^-(x i+ ix l+2 + Xi^iXi^ 2 - 2x i+1 x i - 1 ) . (4) 



The assumption of ergodicity combined with the en- 
ergy constraint Eq. (2) has an important consequence. 
Strikingly, the system exhibits continuous rotational 
symmetry in its statistics, despite the fact that the EOM 
are invariant only under discrete rotations correspond- 
ing to the permutation of indices. The statistics are de- 
termined by only a single parameter, the energy E, or 
cquivalcntly the radius R. 



III. HOPF FUNCTIONAL APPROACH 

Hopf developed a functional method to access the 
equal-time statistics of deterministic dynamical systems. 
In Ref. 2, Hopf derived but did not attempt to solve an 
equation governing the characteristic functional for the 
turbulent velocity field of an incompressible flow. Here 
we briefly review his approach. The reader can find more 
detailed information in Ref. 4. 

For the purposes of illustrating the method, and ap- 
plying it in the next Section to the Orszag-McLaughlin 
system, it suffices to consider deterministic dynamical 
systems described by a set of first-order differential equa- 
tions: 



= Fi(xi,X2, ■ ■ -,x N ); 

at 



1. 



,N. 



(5) 



Introduce a set of variables u conjugate to the coordinates 
x and a functional ^ that resembles a wavefunction in 
quantum mechanics: 



nu,x(t)} 



= e ™-z(*) 



(6) 



It then follows from the EOM, Eq. (5), that is governed 
by a Schrodinger-like equation: 



i — = 
dt 



(7) 



Here the Hopf "Hamiltonian" on the right-hand side is 
the linear operator: 



N 



For odd values of N > 5, Eq. (3) is the only isolating in- 
tegral of motion [5]. In addition, for N = 5, all periodic 
orbits and fixed points are unstable [5]. The system ap- 
pears to be ergodic for all initial conditions except a set 
of measure zero. This also seems to be the case for other 
odd values of N > 5. Finally we note that the dynamical 
system formally obeys Liouville's Theorem, despite being 
non-Hamiltonian: 



Ed ( dxi 
dxl ' ~ 



^ j j \idui' idu 2 '' 

3=1 



1 d 



' i 8un 



(8) 



Note that H is generally non-Hcrmitian, as there is no re- 
quirement (unlike in quantum mechanics) that the time- 
evolution of VP be unitary. For instance the inner prod- 
uct, which plays a central role in quantum mechanics, 
has no particular meaning in the Hopf formalism. We 
call Eq. (7) together with Eq. (8) "Hopf's equation" in 
the following. 

We emphasize that Hopf's equation is linear even 
though the EOM may be nonlinear. Therefore - and this 
is the key point - Eq. (7) can be averaged over an ensem- 
ble of initial conditions. It is the averaging that blurs the 
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deterministic trajectories and, pursuing the analogy with 
quantum mechanics, can be thought of as the origin of 
non-zero H. The corresponding characteristic functional 
t) that solves the Hopf equation is simply an average 
of ty[u, x(t)] over initial conditions at time t = 0: 



'CO 

e lS - s ^P(x(0))d N x(0) . (9) 



For notational simplicity we eliminate the bar over ^ in 
the following. 

It is readily seen that derivatives of the characteristic 
functional yield the equal-time moments: 



1 9* (it, t) 



d 2 ^(u, t) 



dujduk 



x 3 (t)P(x(0))d N x{0) 
Cj(t), (10) 

SO 

x j (t)x k (t)P(S(0))d N x(0) 

OO 

Xj{t)x k (t), (11) 



and so on. All conjugate coordinates u are taken to zero 
after derivatives are taken. From Eq. (9) it is clear that 
the characteristic functional is simply the Fourier trans- 
form of the PDF, and one uniquely determines the other. 
We note that the calculation of n-point correlation func- 
tions are more simply extracted from the characteristic 
functional than from the PDF. In the case of the PDF we 
must integrate over all of the coordinates. For the charac- 
teristic functional all that is required is the calculation of 
some derivatives followed by setting all of the conjugate 
coordinates to zero [6]. The advantage is particularly 
clear in the case of infinite-dimensional problems such 
as fluid flow for which the velocity field is a continuous 
function of position. 

Of special importance are zero-mode characteristic 
functionals that satisfy H^>q = 0. The zero- mode plays 
a role as important as the ground state in quantum me- 
chanics as it encodes information about the steady-state 
behavior of the dynamical system. In the following, we 
employ Dirac's bracket notation (. . .) to indicate steady 
state statistical averages. 

In order for W to be a valid characteristic functional, it 
must obey two general boundary conditions that follow 
directly from Eq. (9): 



*(0) 



1 



and 



lim 

\u\ — >cx 







(12) 



(13) 



Furthermore, a necessary and sufficient condition to en- 
sure semi-positivity of the PDF is specified by Bochncr's 
Theorem [7]: 



d N u d N u' ^{u)^{u-u')(f)*{u') > 



(14) 



for all complex- valued test functions 4>(u). This con- 
strains all of the equal-time moments or cumulants. For 
example, even cumulants cannot be negative [8, 9]. Addi- 
tional boundary conditions arise from conservation laws. 
In particular, for the Orszag-McLaughlin system, con- 
servation of energy imposes an infinite set of boundary 
conditions: 



E r ' 



(15) 



u=0 



for all integer n. 

For completeness, we conclude this brief review of the 
Hopf functional approach by discussing the inclusion of 
random forcing. For simplicity, consider the introduction 
of Gaussian random forcing /j, ^-correlated in time, and 
with zero mean. The EOM are modified to read: 
Ht ■ 

F t (x) + h{t) (16) 



where 



dt 

IfM = o, 

[f {t)h{t')] = 2T6 jk 6(t-t') 



(17) 
(18) 

and the square brackets [. . .] denote an average over the 
random variable. The Hopf equation for the characteris- 
tic functional now reads [10]: 



dt 



u ■ F(-iV u ) - zr|-u|^ f(u, t) 



(19) 



One possible derivation of Eq. (19) begins with the 
better-known Fokker-Planck equation [11]: 



8P(x, t) 
dt 



= [-V x -F(x)+TVi)P{x,t) 



(20) 



The Hopf equation, Eq. (19), is recovered upon making 



the duality transformation i, — > i-J^— and S— _ > I u . 

J J l OUj OXj I ■> 

The duality transformation respects the commutation re- 
lation between the conjugate coordinates: [xj, jUk] = 
8jk- Thus the Hopf and Fokker-Planck approaches are 
seen to be equivalent, and dual to one another, yet the 
inclusion of random forcing seems more intuitive in the 
latter case. As is well-known, the first term in the Fokker- 
Planck equation is a statement of conservation of phase- 
space density. The second is a diffusion term due to ran- 
dom forcing, which smears out an initially sharp distribu- 
tion. Pawula [12] showed that Gaussian random forcing 
leads to only a finite number of terms in the Fokker- 
Planck equation. If the forcing is non-Gaussian, the 
Fokker-Planck equation must include an infinite set of 
terms involving derivatives higher than two. 



IV. EXACT CHARACTERISTIC FUNCTIONAL 

The Hopf Hamiltonian for the Orszag-McLaughlin dy- 
namical system follows directly from Eqs. (1) and (8): 



N 



d 2 



if 



du l+1 du l+2 dui-idiii-2 
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d 2 



du i+ idui 



(21) 



where periodic boundary conditions d/du i+ ^ = d/dui 
are implied. For this particular system the Hamiltonian 
is Hermitian, though (as noted above) this will not always 
be the case for other dynamical systems. 

It is straightforward to show that any ^(u) that is a 
function purely of the magnitude of u is a zero-mode of 
the above H . This is consistent both with the assumption 
of ergodicity, and with the numerical simulation of the 
PDF illustrated in Fig. 1. However, such zero-modes are, 
in general, a superposition of states of different energies, 
Eq. (15). Imposition of the energy constraint Eq. (15) 
picks out the desired solution. Consider a Taylor series 
expansion of m the dimcnsionlcss parameter (i?|u|) 2 : 



(22) 



1=0 



The series coefficients can be found by observing that the 
energy constraint requires 



(23) 



and furthermore, the moments may be expressed in terms 
of the coefficients bi . For example, successive differenti- 
ation of Eq. (22) reveals that 



(xj) = 4! b 2 R\ 



(24) 



and 



{x\xl) = 2 2! 2! b 2 R 4 . (25) 
Combining this with the energy constraint 

i? 4 = N(xt)+2(^y x 2 1 x 2 2 ) 

= NA\ R% + 2 2 ry\ 2! 2! R% (26) 



yields 



iV4! + 2 2 (^)2!2!' 



(27) 



We note that the exact zero-mode characteristic func- 
tional satisfies the boundary conditions, Eqs. (12) and 
(13). 

The characteristic zero-mode may be approached from 
another direction, one that does not rely upon series ex- 
pansion. Knowing that the trajectories are distributed 
uniformly on the hypersphere's surface, the stationary 
PDF is a simple <5-function shell: 



P{x) 



1 



'7V-1 



(30) 



where R is given by Eq. (3). The normalization factor 
l ^ s s i m ply the- inverse of the total surface area of 
the hypersphere. By direct substitution it can be seen 
that this is a zero-mode solution to the Fokker-Planck 
Eq. (20). Now 'J'o may be obtained from the Fourier 
transform of Eq. (30) . Due to the spherical symmetry it is 
most convenient to work in spherical coordinates. The N 
Cartesian coordinates are related to spherical coordinates 
by: 

N-2 

x\ = rsin</)o JJ^ sin^- 
3=1 

N-2 

x 2 = rcosifio Y[ sm <fij 

3=1 
N-2 

£3 = rcos(/>i J^J sin^j ,etc. (31) 

3=2 



In these coordinates the volume element is: 

N N-2 

dv = n dx i = rN ~ idr n shii fa d< t 



(32) 



3=1 



3=0 



By exploiting spherical symmetry, we obtain the charac- 
teristic functional, Eq. (29): 



*o(«) 



1 



e ,M S(\x\ - R) dV 



Sn-i 

T(N/2) (Rlul^) 1 -^ 2 J N/2 ^(R\u\) .(33) 



Because "Jo is the Fourier transform of a non-negative 
PDF, it automatically satisfies Bochner's Theorem, 
Eq. (14). 



Each coefficient can be determined in this manner. The 
series so generated is: 



INFINITE DIMENSIONAL LIMIT 



* (£) = r(iv/2)^ 



e=o 



1 (-1)' 
T(£ + N/2) 2 2 H\ 



(R\u\) 2e . (28) 



The series may be summed into the Besscl function 



J 



N/2-1- 



*o(d) - r(JV/2) {R\u\/2f- N ' 2 J N/2 ^(R\u\) . (29) 



We now show that the statistics are Gaussian in the 
limit of high dimension. Consider the expansion Eq. (28). 
In the limit N — > 00, the Taylor series approaches that 
of a Gaussian: 



00 / 
1=0 



(34) 
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Moments 



Prediction Simulation t\ Simulation t-z 



(xi)/R 
(xl)/R 2 

{x 1 x 2 )/R 2 
(xt)/R 3 

{x\x 2 )/R 3 

(4)/R 4 
(x\xl)/R A 



J_ 





= 0.2 




i 0.0857 
i 0.0286 



0.000610 

0.201 
-0.00165 
0.000505 
-0.000209 
0.0866 
0.0282 



0.0000823 

0.200 
-0.000211 
0.0000945 
-0.0000171 
0.0857 
0.0286 



TABLE I: A comparison of low order moments, scaled by 
R, as calculated from the exact characteristic functional and 
from time averages in a direct numerical simulation. Here 
R = 17.3, ti = 100 and t 2 = 10000. Note convergence of the 
numerical results to the analytical values as the length of the 
time integration increases. 




Alternatively, we may integrate the PDF over all but the 
last Cartesian coordinate: 



P(x N ) 
Since 



R S 



N-2 I 1 _ f XN 
Sn-1 V ^ R 



N-4 



; M < R 



hm ii-f^y 



N-4 



^-(N-4)(x N /R) 2 



(35) 



(36) 



the Gaussian distribution is recovered at large-N. 

The fact that a Gaussian distribution appears in the 
limit N — > oo is not unexpected. Motion along a single 
coordinate is governed by the essentially random values 
of the remaining N — 1 coordinates. The Central Limit 
Theorem then applies, resulting in Gaussian statistics. 



FIG. 2: PDF projected onto a single coordinate. The dots 
are from direct numerical simulation of the N — 5 system; 
solid lines are plots of Eq. (35) for N = 5, 7, and 9, showing 
the approach to Gaussian statistics as N becomes large, and 
the dashed line is a normalized Gaussian with height equal to 
the N=9 PDF. 



Fig. 2 compares this function with the normalized his- 
togram generated from a direct numerical simulation. 
Again there is good agreement between the two meth- 
ods. 



VI. NUMERICAL ANALYSIS 



VII. DISCUSSION 



Finally we make a quantitative comparison with direct 
numerical simulation. To be definite consider the N = 
5 case. The characteristic functional Eq. (29) in this 
specific case reduces to: 



*o(«) = 3 



sin(ii|u|) - R\u\ cos(-R|S|) 



i? 3 | 



or in scries form Eq. (28) to 
" 2 ^| 2 i? 4 | 



*o(«) = 1 - 



R" 



(37) 



(38) 



10 280 15120 

Table I shows that there is good agreement between mo- 
ments calculated by successive differentiation of the exact 
characteristic functional Eq. (37) and the time-averaged 
values obtained from direct numerical simulation. Like- 
wise Eq. (35) gives the projected PDF for the single co- 
ordinate x$: 



P(x B ) 



3 X 5 



x 5 <R 
x 5 > R 



(39) 



We employed the Hopf characteristic functional ap- 
proach to study the statistics of the deterministic Orszag- 
McLaughlin system. The system is simple enough, with 
emergent spherical symmetry in the statistics, that the 
exact characteristic zero-mode can be found. Equal-time 
statistics so obtained are non-Gaussian and are in good 
agreement with direct numerical simulation. 

Exact solutions serve as rigorous illustrations of phys- 
ical principles and as benchmarks against which approx- 
imate solutions can be checked. For instance, in the case 
of the Orszag-McLaughlin system, a cumulant expansion 
carried out through third order reproduces the exact first 
through third cumulants (only the second cumulants are 
non-zero), yet higher-order cumulants are by definition 
zero. Inspection of Eq. (38) immediately reveals the size 
of the errors so incurred. The development of more so- 
phisticated, non-perturbative, methods to address the 
non-equilibrium statistical mechanics of dynamical sys- 
tems can likewise benefit from tests against such exact 
solutions. 
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